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CONVERGENCE RATE FOR A RADAU COLLOCATION METHOD 
APPLIED TO UNCONSTRAINED OPTIMAL CONTROL * 

WILLIAM W. HAGERt, HONGYAN HOLT*, AND ANIL V. RAO § 

Abstract. A local convergence rate is established for an orthogonal collocation method based on 
Radau quadrature applied to an unconstrained optimal control problem. If the continuous problem 
has a sufficiently smooth solution and the Hamiltonian satisfies a strong convexity condition, then 
the discrete problem possesses a local minimizer in a neighborhood of the continuous solution, and 
as the number of collocation points increases, the discrete solution convergences exponentially fast 
in the sup-norm to the continuous solution. An earlier paper analyzes an orthogonal collocation 
method based on Gauss quadrature, where neither end point of the problem domain is a collocation 
point. For the Radau quadrature scheme, one end point is a collocation point. 

Key words. Radau collocation method, convergence rate, optimal control, orthogonal colloca¬ 
tion 


1 . Introduction. A convergence rate is established for an orthogonal collocation 
method applied to an unconstrained control problem of the form 

minimize C(x(l)) 

subject to x(£) = f(x(£), u(£)), £ € [—1,1], (1.1) 

x( 1) = Xo, 

where the state x(£) £ R ra , x = 2 Xj the control u(£) £ R m , f : R" x R m —X R ra , 
C : R” —X R, and x 0 is the initial condition, which we assume is given. Assuming the 
dynamics x(£) = f(x(£), u(£)) is nice enough, we can solve for the state x as a function 
of the control u, and the control problem reduces to an unconstrained minimization 
over u. 

Let V N denote the space of polynomials of degree at most N defined on the 
interval [—1, +1], and let V^ denote the n-fold Cartesian product Vn x ... x Vn- We 
analyze a discrete approximation to m, introduced in ffUdj, of the form 

minimize C(x(l)) 

subject to x(t») = f(x(Tj), Ui), 1 < * < N, (1.2) 

x(-l)=x 0 , 

At the collocation points t», 1 < i < N, the equation should be satisfied. The control 
approximation at time t* is U;. We focus on the Radau quadrature points satisfying 

— 1 < T\ < r 2 < ... < tn = +1. 


The analysis we give also applies to the flipped Radau scheme obtained by reversing 
the sign for each collocation point. In (11.21) the dimension of Vn is N + 1, while there 
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are N + 1 equations in (11.21) corresponding to the collocated dynamics at N points 
and the initial condition. When the discrete dynamics is nice enough, we can solve 
for the discrete state x £ as a function of the discrete controls U;, 1 < i < N, 
and the discrete approximation reduces to an unconstrained minimization over the 
discrete controls. 

In an earlier paper [16] we analyzed a scheme based on Gauss quadrature, where 
the collocation points he in the interior of the interval [—1,+1], The Gauss scheme 
is easier to analyze than the Radau scheme of this paper due to the symmetry of the 
Gauss collocation points, and the fact that the none of the collocation points lies at 
an end of the problem domain. For the Radau scheme, symmetry is broken and the 
polynomials in the discrete adjoint equation have degree N— 1 compared to the degree 
N polynomials used for the state approximation. As will be seen, the presence of the 
Radau collocation point at the end of the interval [—1,+1] leads to the embedding 
of the terminal adjoint condition of the discrete problem into the discrete adjoint 
dynamics. Despite these differences in the analysis and despite the fact that Gauss 
quadrature has a higher degree of accuracy than Radau quadrature, the convergence 
rate obtained for the Radau scheme is exactly the same as that of the Gauss scheme. 
Moreover, in numerical experiments with test problems where an exact solution is 
known, the observed error in the state, control, and adjoint for the Radau scheme 
is very similar to the observed error for the Gauss scheme. The fact that a Radau 
quadrature point can be placed at the end point of the interval leads to a simpler 
implementation of terminal constraints and terminal cost. And the Radau scheme is 
easier than the Gauss scheme to extend to an /ip-framework [11111201111] where the 
interval [—1, +1] is partitioned into a mesh and a different polynomial is employed in 
each mesh interval. 

The analysis in this and the earlier paper m needs further extensions in order to 
handle Lobatto collocation schemes such as those in HIS] where iq = —1 and tjv = +1. 
Although the Gauss and Radau scheme have similar errors, the Lobatto scheme can 
converge at a slower rate, as observed in [T2]. An advantage of the Lobatto scheme is 
that the value of the optimal control is estimated at the initial point t = — 1 and the 
terminal point t = +1. Moreover, both initial and terminal constraints are easier to 
implement in the Lobatto framework. A consistency result for a Lobatto collocation 
scheme applied to optimal control is given in m- Other quadrature points that have 
been exploited in the optimal control literature include the Chebyshev quadrature 
points mm, and the extrema of Jacobi polynomials [26 1. 

Our goal is to show that if (x*, u*) is a local minimizer for O), then the discrete 
problem (11.21) has a local minimizer (x jV , u N ) that converges exponentially fast in N 
to (x*, u*) at the collocation points. Convergence rates have been obtained previously 
when the approximating space consists of piecewise polynomials as in m nun m nn 
QH23]. In these earlier results, convergence is achieved by letting the mesh spacing 
tend to zero. In our results, on the other hand, convergence is achieved by letting N, 
the degree of the approximating polynomials, tend to infinity. 

Let C fc (M") denote the space of k times continuously differentiable functions x : 
[—1, +1] —> R ra with the sup-norm || • Hoc given by 

Moo = sup{|x(f)| : t £ [-1, +1]}, (1.3) 

where | -1 is the Euclidean norm. It is assumed that (11.11) has a local minimizer (x* , u*) 
in C 1 (M n ) x C 0 (R m ). Given y £ R n , the ball with center y and radius p is denoted 

&p{y) = {x £ K" : |x — y| < p\. 
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It is assumed that there exists an open set f2 C R” l+n and p > 0 such that 
B p (x*(t), u*(£)) C O for all t £ [—1,+1]. 

Moreover, the first two derivative of / and C are continuous on the closure of fl and 
on £> p (x*(l)) respectively. 

Let A* denote the solution of the linear costate equation 

A*(t) = -V aJ jy(x*(t),u*(t),A*(t)), A*(l) = VC(x*(l)), (1.4) 

where H is the Hamiltonian defined by H(x, u, A) = A T f(x, u). Here VC denotes 
the gradient of C. By the first-order optimality conditions (Pontryagin’s minimum 
principle), we have 


V„ff(x*(i),u*(i),A*(i)) = 0 


(1.5) 


for all t £ [—1, +1]. 

Since the discrete collocation problem (11.21) is finite dimensional, the first-order 
optimality conditions (Karush-Kuhn-Tucker conditions) imply that when a constraint 
qualification holds [2Tj . the gradient of the Lagrangian vanishes. By the analysis in 
mna, the gradient of the Lagrangian vanishes if and only if there exists A £ V^ l _ 1 
such that 


A 0 = A(—1) (1.6) 

A(Tj) = -V x H‘(x(r i ),u i , A(rj)), 1 < i < N, (1.7) 

A(l) = -V X H (x(l), UjV) A(l)) + (A(l) - VC(x(l))K, (1.8) 

0 = X/ U H (x(r,), Ui, A(ri)), 1 <i< N, (1.9) 


where u>i is the Radau quadrature weight associated with Ti , and Ao is the multiplier 
associated with the initial condition x(—1) = Xo in (11.21) . Note that in m, CD is 
written in the form 


N 

VC(x(l)) = A 0 - ^WiVzU(x(Ti),u,, A(r,)). 

i—1 

However, utilizing (EZD. (HUD, and the fundamental theorem of calculus, this reduces 
to the more compact form m - 

In comparing the first-order conditions for Radau collocation to the first-order 
conditions for Gauss collocation nn, the differences are that in Gauss collocation, 
A £ not V^-i- Moreover, in Gauss collocation, the terminal condition for the 
discrete adjoint is simply A(l) = VC(x(l)), while in Radau collocation, the discrete 
costate dynamics and the terminal condition are mixed together as in (11.81) . 

The assumptions utilized in the convergence analysis are the following: 

(Al) x* and A* £ C vJrl for some rj > 3. 

(A2) For some a > 0, the smallest eigenvalue of the Hessian matrices 
V 2 C(x*(l)) and V^ u) i7(x*(f),u*(f),A*(f)) 
is greater than a , uniformly for t £ [—1, +1]. 
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(A3) The Jacobian of the dynamics satisfies 

||V a: f(x*(t),u*(i))|| 00 < 1/4 and ||V I f(x*(i),u*(i)) T || 00 < 1/4 

for all t £ [—1, +1] where || • ||oo is the matrix sup-norm (largest absolute row 
sum), and the Jacobian V^f is an n by n matrix whose *-th row is (Vj/i) T . 

The smoothness assumption (Al) is used to obtain a bound for the accuracy 
with which the interpolant of the continuous state x* satisfies the discrete dynamics. 
The coercivity assumption (A2) ensures that the solution of the discrete problem 
is a local minimizer. The condition (A3) enters into the analysis of stability for 
the perturbed dynamics; this condition can be eliminated when the problem domain 
[—1, +1] is partitioned into subintervals with a different polynomial on each subinterval 
For the global polynomials analyzed in this paper, (A3) could 
be replaced by any condition that ensures stability of the linearized dynamics. 

In addition to the 3 assumptions, the analysis utilizes 4 properties of the Radau 
collocation scheme. Let tq = — 1, a noncollocated point, and define 


N 

Dij = 1 < i < N, 0 < j < N, where Lj{r) := TT - (1.10) 

f = p T 3 - A 

Here the dot denotes differentiation, and D is a differentiation matrix in the sense 
that (Dp); = p(Ti), 1 < * < N, whenever p £ Vn is the polynomial that satisfies 
p(Tj) = pj for 0 < j < N. The submatrix Di,jv consisting of the tailing N columns 
of D has the following properties: 

(PI) D x : jv is invertible and ||D 1 .] V || 00 = 2. 

(P2) If W is the diagonal matrix containing the Radau quadrature weights u; on 
the diagonal, then the rows of the matrix [W 1 ' 2 Di : at] _ 1 have Euclidean norm 
bounded by \/2. 

The fact that Di : jv is invertible is established in m Prop. 1], and a formula for the 
elements of is given in m equation (53)]. From the formula, the elements in the 
last row of Djf N are the Radau quadrature weights, which are positive and sum to 2. 
Although elements in the earlier rows of can be either positive or negative, we 

find numerically that their absolute sum is always less than 2. Similarly, the elements 
in the last row of [W^Dnjv] -1 are the square roots of the Radau quadrature weights. 
Hence, the Euclidean norm of the last row of [W 1,/2 Di : ]v] _ 1 is \/2. Numerically, we 
find that Euclidean norm of the earlier rows is always less than \/2. 

Let by the N by N matrix defined by 


D h = 



1 < i < N, 


1 <j<N. 


The matrix D* arises in the analysis of the costate equation. In Section 4.2.1 of [T2] , 
we introduce a matrix which is a differentiation matrix for the collocation points 
Ti, 1 < i < N. That is, if p is a polynomial of degree at most N — 1 and p is the 
vector with components p(ri), 1 < i < N, then (D^p)j = p(ri). The matrix D* only 
differs from in a single entry: D^ NN = D^ NN — 1 /u>n- As a result, 

(D t p) i = p(Ti), 1 < i < N, (D*p)jv = p(t n ) - p{l)/oj N . (1.11) 

If Dip = 0, then j)(r,;) = 0 for i < N by (11.111) . Since p has degree N — 2 and it 
vanishes at iV— 1 points, p is identically zero and p is constant. By the final equation 
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in (HUD, p(l) = 0 when D^p = 0, which implies that p is identically zero. This shows 
that D* is invertible. We find that D* has the following properties: 

(P3) D* is invertible and ||(D *) _1 || 00 < 2. 

(P4) The rows of the matrix [W^D*] -1 have Euclidean norm bounded by 
In Proposition 19.II at the end of the paper, an explicit formula is given for the inverse 
of D*. However, it is not clear from the formula that ||(Dt) 1 || 00 is bounded by 
2. Numerically, we find that the norms in (P3) and (P4) achieve their maximum in 
the first row of the matrix, and these norms increase monotonically towards the given 
bounds. Properties (P1)-(P4) differ from the assumptions (Al)-(A3) in the sense that 
the properties seem to hold for any choice of N, although a proof is missing, while 
(A1)-(A3) only hold for certain control problems. In the analysis of the Gauss scheme 
m, properties (P3) and (P4) followed immediately from (PI) and (P2) since it could 
be shown that the discrete costate matrix was related to the state differentiation 
matrix through an exchange operation. However, due to the asymmetry of the Radau 
collocation points and the lower degree of the polynomials in the discrete adjoint 
system disD dp , the relation between the state and costate matrices for the Radau 
scheme is not clear. Nonetheless, the bounds in (P3) and (P4) are observed to be the 
same as the bounds in (PI) and (P2). 

If x N is a solution of (11.21) associated with the discrete controls Uj, 1 < * < IV, 
and if X N £ satisfies (11.61) (11.91) . then we define 


X N =[ 

x "( r o)> 

x "(ti), ..., 

X W (tat) 

X* = l 

x *(t 0 ), 

x*(ti), ..., 

X*(tat) 

U N =[ 


Ui, 

u N 

U* = [ 


u*(n), ..., 

u*(t/v) 

A n =[ 

A" (t 0 ), 

A "(n), ..., 

A n (t n ) 

A* = [ 

A* (t 0 ), 

A*(ti), ..., 

A* (tat) 


For any of the discrete variables, we define a discrete sup-norm analogous to the 
continuous sup-norm in (11.31) . For example, if U N £ R mJV with U i £ R m , then 

IIU^Hoo = sup{|Uj| : 1 <i < N}. 

The following convergence result is established: 

Theorem 1.1. //(x*,u*) is a local minimizer for the continuous problem (ED 
and both (A1)-(A3) and (P1)-(P4) hold, then for N sufficiently large with N > 77 -1-1, 
the discrete problem ED has a local minimizer (X N ,U N ) for which 

max {||X" - Xl^ HU" - U’lU, ||A" - A*|U} < cN 2 ~\ (1.12) 

where c is independent of N. 

The discrete problem provides an estimate for optimal control at r/v = +1, how¬ 
ever, there is no discrete control at ro = — 1. Due to the strong convexity assumption 
(A2), an estimate for the discrete control at t = — 1 can be obtained from the min¬ 
imum principle ED since we have estimates for the discrete state and costate at 
tq = — 1. Alternatively, polynomial interpolation could be used to obtain estimates 
for the optimal control at t = — 1. 

The paper is organized as follows. In Section [2] the discrete optimization problem 
ED is reformulated as a nonlinear system of equations obtained from the first-order 
optimality conditions, and a general approach to convergence analysis is presented. 
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Section [3] obtains an estimate for how closely the solution to the continuous prob¬ 
lem satisfies the first-order optimality conditions for the discrete problem. Section [J] 
proves that the linearization of the discrete control problem around a solution of the 
continuous problem is invertible. Section [5] establishes an L 2 stability property for 
the linearization, while Section [6] strengthens the norm to L°°. This stability prop¬ 
erty is the basis for the proof of Theorem 11.11 A numerical example illustrating the 
exponential convergence result is given in Section 0 

Notation. The meaning of the norm || • Hoc is based on context. If x € C°(R"), 
then Hxlloo denotes the maximum of |x(t)| over t € [—1, +1], where | • | is the Euclidean 
norm. If A £ R mx ", then HAHoo is the largest absolute row sum (the matrix norm 
induces by the vector norm). If U £ M. mN is the discrete control with Uj £ R m , 
then IIUIU is the maximum of |Uj|, 1 < i < N. The dimension of the identity matrix 
I is often clear from context; when necessary, the dimension of I is specified by a 
subscript. For example, I„ is the n by n identity matrix. VC denotes the gradient, 
a column vector, while V 2 C denotes the Hessian matrix. Throughout the paper, c 
denotes a generic constant which has different values in different equations. The value 
of this constant is always independent of N , the degree of the polynomials used in the 
discrete approximation of the state. 1 denotes a vector whose entries are all equal to 
one, while 0 is a vector whose entries are all equal to zero, their dimension should be 
clear from context. If D is the differentiation matrix introduced in (11.101) . the Dj is 
the j-th column of D and 'Dtj is the submatrix formed by columns i through j. 

2. Abstract setting. Given x £ V JJ- and u £ R mN that are feasible in (11.21) . 
define X, = x(rj) and Uj = Uj. As shown in [T2], the discrete problem (11.21) can be 
reformulated as the nonlinear programming problem 

minimize C(Xj\r) 

subject to Ef=o Aj X , = ffX.A,), 1 <i < N, (2.1) 

Xo = xo. 

Also, m shows that the equations obtained by setting the gradient of the Lagrangian 
to zero are equivalent to the system of equations 

N 

A 0 = VC(Xv) + w i V x fl r (X i , U i; Aj), (2.2) 

2—1 
N 

A j = -V x H{X i ,U i ,A i ) t 1 <i<N, (2.3) 

3 =1 

N 

J2 = -V X H (X N , Ujv, A n ) - X7 C(X n )/uj n , (2.4) 

j=i 

0 = V u fr(X i ,U i ,A < ), 1 <i<N, (2.5) 

where Ao is the multiplier associated with the equation Xo = Xo and Aj for i > 0 is 
related to the Lagrange multiplier Aj associated with the *-th equation in the discrete 
dynamics by 

Aj = Aj/wj. (2.6) 

The first-order optimality conditions for the nonlinear program m consist of 
the equations (12.21) (12.51) , and the constraints in (12.11) . This system can be written as 
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T(X, U, A) = 0 where 

(71,72, ■ • •, 7fs)(X, U, A) e R niV xl"xl"x ^n(iv-i) x R n x 
The 6 components of T are defined as follows: 

Tii(X, U, A) = ^ AjXjj - f(Xj, Uj), 1 < z < AT, 

T 2 (X,U,A) = X 0 -xo, 

N 

r 3 (X,U, A) = A 0 - VC(Xjv) - Ai), 

r 4l (X,U,A) = +V x H(X i ,V i ,A i ), I < i < N, 

N 

r 5 (X, U, A) = D%j A j + (Xjv, Ujv, A n ) + X7C(X n )/lu n , 
i=i 

r 6 ,;(X,U, A) = V„H(X i ,U i) AO, l<i<N. 

The proof of Theorem 11. 1 1 reduces to a study of solutions to T(X, U, A) = 0 in 
a neighborhood of (X*, U*, A*). Our analysis is based on 0 Proposition 3.1], which 
we simplify below to take into account the structure of our T. Other results like this 
are contained in Theorem 3.1 of [Tj, in Proposition 5.1 of El. and in Theorem 2.1 of 

El- 

Proposition 2.1. Let X be a Banach space and y be a linear normed space with 
the norms in both spaces denoted || • |. Let T: X \—> y with T continuously Frechet 
differentiable in B r (9*) for some 9* £ X and r > 0. Suppose that 

|| VT(0) - VT(0*)|| < e for all 9 G B r (0*) 

and VT(9*) is invertible; and define p := ||VT(#*) -1 ||. If sp < 1 and \\T (0*)|| < 
(1 — pe)r/p, then there exists a unique 0 G B r (9*) such that T(9) = 0. Moreover, we 
have the estimate 

||0 - 0*|| < — ||T (0*)|| ■ (2.7) 

1 — pe 

We apply Proposition O with 9* = (X*,U*,A*) and 9 = {X N ,U N ,A N ). The 
key steps in the analysis are the estimation of the residual ||T(0*)||, the proof that 
VT(0*) is invertible, and the derivation of a bound for ||VT(0*) -1 || that is indepen¬ 
dent of N. In our context, the norm on X is 

||0|| = ||(X,U, A)|U = maxdIXIU, HUH*,, UA^}. (2.8) 

For this norm, the left side of (11.121) and the left side of (12.71) are the same. The 
norm on y enters into the estimation of both the residual ||T(0*)|| in (12.71) and the 
parameter /i:=||VT(0*) -1 ||. In our context, we think of an element of y as a vector 
with components y», 1 < i < 3 N + 2, where y* G M n for 1 < i < 2 N + 2 and y» G 
for i > 2 N + 2. For example, 7i(X,U, A) G W lN corresponds to the components 
yi G R n , 1 < i < N. For the norm in y, we take 

IIyIIoo = sup{|y,| : 1 < i < 3N + 2}. 


(2.9) 
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3. Analysis of the residual. We now establish a bound for T(X*, U*, A*), the 
residual which appears on the right side of the error bound (12.71) . This bound for the 
residual ultimately appears in the right side of the error estimate (11.121) . 

Lemma 3.1. //(Al) holds, then there exits a constant c, independent of N, such 

that 

HT^inA^IU^cA 2 -” (3.1) 


for all A > rj + 1. 

Proof. By the definition of T, we have 

r 2 (X*, u*, A*) = X* - x 0 = x*(r 0 ) - x 0 = x*(-l) - X 0 = 0 

since x* satisfies the initial condition in (11.11) . Likewise, 7(3(X*, U*, A*) = 0 since 
m holds for all t £ [— 1,+1], which implies that (11.51) holds at the collocation 
points. 

Now consider 71 ■ Since D is a differentiation matrix associated with the colloca¬ 
tion points, we have 

N 

J2 A;. 7 X* = x j (tO, 1 < * < A, (3.2) 

j=o 

where x J £ is the (interpolating) polynomial that passes through x* (jj) for 0 < 
j < A. Since x* satisfies the dynamics of m, 

f(X*, U*) = f(x*(Ti), u*(t,)) = x*(n). (3.3) 

Combine (13.21) and (13.31) to obtain 

Tu(X*,U*,A*) =x I (T i )-x*(r i ). (3.4) 

Proposition 2.1 and Lemma 2.2 in m yield 

||i 7 - i'lloo < [(1 + 2A 2 ) + 6eA(l + Cl log N)] 

for all A >77 + 1 , where x^ +1 ) is the ( 77 + l)-st derivative of x and Ci log A is a bound 
for the Lebesgue constant of the point set Tj, 0 < j < A, given in Theorem 2.1 of 
[25] . Hence, there exists a constant c 2 , independent of A but dependent on 77 , such 
that 

|Tii(X*,U*,A*)| = |x J (Ti) -x*(rO| < ||x 7 -x^U < c 2 N 2 ~\ (3.5) 

which complies with the bound (13.11) . 

Next, let us consider 71- By (11.111) if X 1 £ is the (interpolating) polynomial 

that passes through A *( tj ) for 1 < j < A, we have 

N N 

A^A^rO. 1 <*<A, J2 Dt Nj a; = A J (1) - A J (1)/ W 3V. (3.6) 

1 =1 i =1 

Since A* satisfies (ED, it follows that 


V x A(X*,U*,A*) = V x 77(x*(r i ),u*(r i ),A*(r i )) 


a*(a). 


(3.7) 
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Hence, for i < N, we have 

7i(X*,U*,A*) = A 7 (t0 - A *(Ti). (3.8) 

By a similar analysis as that used for 71 in (13.51) . we conclude that 

|T 4 i(X*,U*, A*)| < cN 2 - 1 ! for i < N. (3.9) 

The difference between the analysis of the state in (13.51) and the analysis of the costate 
in m is that the 0(\ogN) bound for the Lebesgue constant of r 0 , ..., must 
be replaced by an 0(y/N) bound, derived in Theorem 5.1 of El, for the Lebesgue 
constant of ti, ..., tjv- This difference between the state and the costate arises 
since the state interpolant x 7 G interpolates x* at To, ..., tjv while the costate 
interpolant A 7 G 7 ? ](r_ 1 interpolates A* at n, ..., tn■ Since the Lebesgue constant 
term in the bound for ||71i(X*, U*, A*)|| is dominated by the other terms, the right 
sides of m and m are the same. 

Similarly, using (13.71) and the second equation in (13.61) . it follows that 

r 5 (X*,U*, A*) = A 7 (l) - A*(l) + (VC(x*(l)) - A*(1))/ WJV = A 7 (l) - A*(l) 

since A 7 (l) = A*(l) = VC(x*(l)) by (11.41) . Hence, just like the bound for 71,; in (13.91) . 
we have 

|7' 5 (X*,U*,A*)| < cN 2 ~ v . 

Now consider 71- By (13.71) . the definition Ag = A*(—1), and the terminal condi¬ 
tion A*(l) = VC(x*(l)) from (11.41) . we have 

N 

75(X*,U*,A*) = A*(-l)-A*(l) + 53wiA*(Ti). (3.10) 

2 — 1 

By the fundamental theorem of calculus and the fact that IV-point Radau quadrature 
is exact for polynomials of degree up to 2 N — 2, we have 

r i n 

0 = A 7 (-l) - A 7 (l) + / X 1 (t)dt = A 7 (—1) — A 7 (l) + ujjX 1 (jj). (3.11) 

J - 1 i=i 

Subtract (|3.11f) from ()3.10[) to obtain 

N 

71(X*, U*, A*) = A*(—1) - A 7 (-l) + £>,- (a* (tj ) - A 7 (tj)) (3.12) 

j =i 

since A 7 (l) = A*(l). Since Wj > 0 and their sum is 2, it follows from (13.81) and (13.91) 
that 

N 

U}j 

3 = 1 

By Theorem 15.1 in |24j and Lemma 2.2 and Theorem 5.1 in E, we have 

IA*(—1) - A 7 (—1)| < (1 + Cl VN) " +1 1| A*^ |U 

< dV- (0 ' 5+7?) . (3.14) 



A'fo)- X*( Tj ) 


< cN 2 ~ v . 


(3.13) 


We combine (13.121) (13.141) to see that 7s(X*,U*,A*) also complies with the bound 
(EH). This completes the proof. 0 
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4. Invertibility. In this section, we show that the derivative VT(0*) is invert¬ 
ible. This is equivalent to showing that for each y £ y, there is a unique 9 £ X 
such that VT(0*)[6»] = y. In our case, 9* = (X*,U*,A*) and 9 = (X,U,A). To 
simplify the notation, we let VT*[X, U, A] denote the derivative of T evaluated at 
(X*, U*, A*) operating on (X, U, A). This derivative involves the following 6 matri¬ 
ces: 

Aj = V x f(x*(Tj), U*(Tj)), Bj = V„f(x*(7i),U*(Tj)), 

Q; = (x*(r»), u*(r 4 ), A*(r»)), S» = S7 XU H (x*(t»), u*(Tj), A*(rj)), 

R-i = ^uuH (x*(r»), u*(r»), A*(ri)), T = V 2 C'(x*(l)). 

With this notation, the 6 components of VT*[X, U, A] are as follows: 

VTh[X,U,A] 

vt;[x,u,a] 
vt 3 *[x,u,a] 

vr 4 *[x, u,a] 

vt 5 *[x,u,a] 
vt 6 *[x,u,a] 

The analysis of invertibility starts with results concerning the invertibility of the 
linearized discrete state dynamics. 

Lemma 4.1. If (PI) and (A3) hold, then for each q £ K" and p £ K nAr with 
Pi £ R n , the linear system 

N \ 

Y DijXj - A,X, = pi 1 < i < N, (4.1) 

j=o ) 

X 0 = q, (4.2) 

has a unique solution Xj £ K", 0 < j < TV. This solution has the bound 

11 Xy 11 oo < 4||p|| oo + 2|| q|| oo, 0 < j < N. (4.3) 

Proof. If X is a solution of (14.11) (14.21) . then Xo = q and Xo trivially satisfies 
(TOl) . Next, focus on the remaining components of X. Let X be the vector obtained 
by vertically stacking X 3 through Xjv, let A be the block diagonal matrix with i-th 
diagonal block A.;, 1 < i < N, and define D = Di : jv <8> I n and Do = Do ® I„ where <g) 
is the Kronecker product. With this notation, the linear system (roi dxa reduces 
to 


— ^ 'y ' — AiXj BiUi, 1 < * < A, 

= X 0 , 

N 

= A 0 - TXjv - ]Twi(A7A ?; + Q*Xi + SiUO, 

i=1 

= D h A ^j + A 7 A * + QiXi + SlUl , 1 <i<N, 

= ( ^A'-j A/^j + A](rAiv + QnXjv + SivUiv + TXjv/waT) 

= SjX, + R,U, + BjA x , 1 < i < N. 


(D A)X = p D 0 q. 


(4.4) 
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By (PI), Dx ; jv is invertible which implies that D is invertible and D 1 = D-^ ® In- 
Since the polynomial that is identically equal to 1 has derivative 0 and since D is a 
differentiation matrix, we have D1 = 0, which implies that D^Do = —1. It follows 
that 


D -1 Do = [D-jy. ® I„][D 0 8) I„] = -1 ® In- 
Multiply (Id.41) by D _1 to obtain 


(I — D _1 A)X = D -1 p + (1 ® I n )q. (4.5) 

By (PI) HD^jylloo^ 2, which implies that ||£)~^||oo < 2. By (A3) HAH^ < 1/4. By 
[151 p. 351], I — D _1 A is invertible and ||(I — D _1 A) _1 || oo < 2. Multiply (14.51) by 
(I — D -1 A) -1 and take the norm of each side to obtain HXH^ < dHpH^ + 2||q||oo. 
This complete the proof of (14.31) . □ 

Next, we establish the invertibility of VT*. 

Proposition 4.2. If (PI), (A2), and (A3) hold, then VT* is invertible. 

Proof. Our approach is to formulate a strongly convex quadratic programming 
problem which has a unique solution (X, U) by (A2), and which has the property 
that the associated first-order optimality condition is VT*[X,U, A] = y. Since VT* 
is square and VT*[X, U, A] = y has a solution for each choice of y, we conclude the 
VT* is invertible. 

The quadratic program is 


minimize ±Q(X, U) + £(X, U) 

subject to V/ 0 /h,X, - A,X, - 15,U, + y,,. 1 < i < N, ’ 

Xo = Y2, 

where the quadratic and linear terms in the objective are 

N 

Q(x,u) = x J N tx n + ( x 7 Q* x * + 2X7S.U,; + u^ir), 

2—1 

( N \ N 

Y3+J2 - E (y T ^+■ 

2=1 J 2=1 


(4.6) 

(4.7) 

(4.8) 


In (14.8[) we simplified the formula for C by introducing y 4 ^v = ys- By Lemma IXT1 the 
quadratic programming problem (14.61) is feasible. Since the Radau quadrature weights 
u>i are strictly positive, it follows from (A2) that Q is strongly convex. Hence, there 
exists a unique optimal solution to (14.61) for any choice of y. Since the constraints are 
linear, the first-order optimality conditions hold. The linear term was chosen so that 
the first-order optimality conditions for (14.61) reduce to VT*[X, U, A] = y. Hence, 
the existence of a solution to (14.61) for each y would imply that VT*[X, U, A] = y 
has a solution for each choice of y. To complete the proof, we need to show that the 
first-order optimality conditions for (14.61) are equivalent to VT*[X,U, A] = y. 

The Lagrangian of (14.61) is given by 


N 


- Q(X, U) + £(X, U) + ^ A/ A,X 4 + BiUi 


N 

■ yu - ^2 D 

3=0 




+ Aj(y 2 -X 0 ). 
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The first-order optimality conditions are obtained by setting to zero the derivative 
of the Lagrangian with respect to each of the components of X and U. We give the 
derivation of the 3rd, 4th, and 5th components of VT*[X, U, A] = y, while the 6th 
component follows in a similar fashion, and the 1st and 2nd components are simply 
the constraints in (14.61) . 

Setting to zero the partial derivative of the Lagrangian with respect to X*, 1 < 
i < N, yields the equation 


N 

A i \i T ujjQ^X^ T cjjSjUj ^ ^ Dji\j — lu ^y 4 2■ 

j =i 

Substituting Dji = — ( uii/Uj)Dj '■ and A_, = uijAj, we obtain 

^E D lj A ^ + A.JAi + Q,X, + SiUi = y 4i , (4.9) 

which gives the 4th component of VT* [X, U, A] = y. In a similar fashion, setting 
to zero the partial derivative of the Lagrangian with respect to Xjv yields the 5th 
component of VT*[X, U, A] = y: 

(E -Pjr,- Aj^ + A^Aat + QjvXv + SatUat + TX^v/wat = y 4 N- (4-10) 

Setting to zero the partial derivative of the Lagrangian with respect to Xo gives 
the equation 


N 

A 0 + E (AjDjo - ujjyn) = y 3 . (4.11) 

*=l 

Since D is a differentiation matrix and D1 = 0, it follows that 

N 

DiO — ^ ) Hf j ■ 

1=1 


Consequently, we have 

N N N N N N N 

E^oa, = - EE^Ea = EE w 4 A r 

i— 1 i= 1 j= 1 2=1 j—1 2=1 j=l 

We make this substitution as well as (14.91) and (14.101) into (14.111) . The terms cancel 
to give 

N 

Ao — TXjv — E w *( j E-A* + QjXj + SiUO = y3, 
which is the 3rd component of VT*[X, U, A] = y. This completes the proof. □ 
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5. w-norm bounds for the state and control. In this section we obtain a 
bound for the (X, U) component of the solution to VT*[X, U, A] = y in terms of 
y. Since Xo must satisfy the constraint Xo = y 2 , it is trivially bounded in terms of 
||ylloo- Hence, we focus on (Xj,Uj), 1 < i < N. The bound we derive in this section 
is in terms of the w-norms defined by 

N N 

||X||2 =|X^| 2 + ^ Wl |X i | 2 and ||U|| 2 = ^^|U,| 2 . (5.1) 

2=1 2=1 

This defines a norm since the Radau quadrature weight coi > 0 for each i. Since the 
(X, U) component of the solution to VT*[X, U, A] = y is a solution of the quadratic 
program (HU), we will bound the solution to the quadratic program. 

First, let us think more abstractly. Let ir be a symmetric, continuous bilinear 
functional defined on a Hilbert space T-L , let £ be a continuous linear functional, let 
<j> £ H, and consider the quadratic program 

min | -ir(v + </>, v + (/)) + £{v + (j)) : v £ V j , 

where V is a subspace of 'H. If w is a minimizer, then by the first-order optimality 
conditions, we have 

n(w, v ) + v ) + £{v) = 0 for all v £ V. 

Inserting v = w yields 


tv(w,w) = —( 7 t(w,4>) +£(w)). 


(5.2) 


We apply this observation to the quadratic program (Id.61) where we treat X 0 = y 2 
as fixed, so the minimization is over (Xj,Uj), 1 < i < N. We identify £ with the 
linear functional £ in (|4.8D but with the Xo term dropped since it is fixed, and 7 r with 
the bilinear form associated with the quadratic term (ED- The subspace V is the 
null space of the linear operator in (Id.61) and (f> is a particular solution of the linear 
system. The complete solution of (14.61) is the particular solution plus the minimizer 
over the null space. 

In more detail, let x denote the solution to (ED (S2I given by Lemma 14.11 for 
p = yi and q = y 2 - We consider the particular solution (X, U) of the linear system in 
(14.61) given by (x,0). The relation (15.21) describing the null space component (X, U) 
of the solution is 


Q(x,u) 


XaTX 


N ' 


N 

E 


[(Qi%i y4i) Xj 



(5.3) 


Here the terms containing x are associated with 7 r(w,<j)), while the remaining terms 
are associated with £, or equivalently with £. By (A2) we have the lower bound 

2(X, U) > a(||X|| 2 + ||U|| 2 ). (5.4) 


All the terms on the right side of (| 5.3 1) can be bounded with the Schwarz inequality; 
for example, 

N / N \ 1 / 2 / N 

E wi xL x * < (E^i^i 2 ) ( E^i x *i 2 

i= 1 \i=l J \i=l 

< V2||y 4 ||oo(||X|| 2 +||U|| 2 ) 1/2 . 


(5.5) 
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The last inequality exploits the fact that the oji sum to 2 and |y 4 » | < ||y 4 ||oo- To 
handle the terms involving x i n (15.31) . we utilize the upper bound Hxjlloo < bUyll^ 
based on Lemma 14.11 with p = yi and q = y 2 ■ Combining upper bounds of the form 
(15.51) with the lower bound (15.41) . we conclude from (15.31) that both ||X|| U and ||U|| W 
are bounded by a constant times ||y||oo- The complete solution of (14.61) is the null 
space component that we just estimated plus the particular solution (x, 0). Again, 
since ||xj||oo < 6||y ||oo > we obtain the following result. 

Lemma 5.1. If (A2)-(A3) and (PI) hold, then there exists a constant c, inde¬ 
pendent of N, such that the solution (X, U) of (14.61) satisfies HXH^ < c||y||oo and 

IIUL < c||y|U- 

6. oo-norm bounds. We now need to convert these w-norm bounds for X and 
U into oo-norm bounds and at the same time, obtain an oo-norm bound for A. As 
in Lemma im the solution to the dynamics in (14.61) can be expressed 

X = (I - D -1 A) -1 [D _1 BU + D-'yi + (1 ® I„)y 2 ] , (6.1) 

where B is the block diagonal matrix with i-th diagonal block B^ and U is obtained 
by vertically stacking Ui through U n- By Lemma FTT1 

||(I - D^ 1 A) _1 [D _1 yi + (1 ® I„)y 2 ] < AWy^ + 2^21^. (6.2) 

The term D : BU can be bounded using (P2) and the strategy given in Section 6 of 
m- That is, we first observe that 

D _1 BU = [D-^ ® I„]BU = [(W 1//2 D 1: jv) _ 1 ® InjBU,, (6.3) 

where W is the diagonal matrix with the quadrature weights on the diagonal and is 
the vector whose i-th element is y/oJi LJ*; the y/UJi factors in (16.31) cancel each other. An 
element of the vector D _1 BU is the dot product between a row of (W 1 / 2 Di : jv) _1 ®In 
and the column vector BU U . By (P2) the rows of (W 1 / 2 Di : at) _ 1 ® I n have Euclidean 
length bounded by y/2. By the properties of matrix norms induced by vector norms, 
we have 


||BU„|| 2 < ||B|| 2 ||U„|| 2 = ||B|| 2 ||U|U. 

It follows that 

||D- 1 BU|| 00 < v / 2||B|| 2 ||U|| LJ < c||yHoc, 

where the generic constant c is independent of N by Lemma l5.ll Combining this with 
( 16 . 21 ) , we conclude that for some constant c, the X component of the solution to ( 14 . 61 ) 
satisfies HXH^ < c||y||oo. Since |X 0 | = |y 2 1 < ||y||oo, the entire X-component of the 
solution to ( 14 . 61 ) satisfies HXH^ < c||y||oo for some c. 

Next, we focus on the 4th and 5th components of VT*[X, U, A] = y which can 
be written 


D*A + A t A + QX + SU + { —) (e N ® I„)TXjv = y 4 , (6.4) 

\UN ) 

where <%) I n , A is obtained by vertically stacking Ai through An, Q and 

S are block diagonal matrices with i-th diagonal blocks Q; and S; respectively, and 
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ejv S R w is the vector whose components are all zero except for the A-th component 
which is 1. Similar to our manipulations of the state dynamics in (16.11) . we use (A3) 
and (P3) to solve for A: 


A = -(I + D t ' 1 A T )" 1 D t - 1 


SU + QX + 



(e^v ® I n )TXjv — y4 


(6.5) 


Let p be the polynomial that is identically one, and let p be the associated 
vector whose components are all one. Making this substitution in (11.111) gives D+l = 
—ejv/w/v, which implies that 

D* _1 ejv = — watI. (6.6) 


Consequently, we have 

fi* _1 (ejv ® I n )/wN = [D* _1 ® I n ](ejv ® I„)/wjv = —1 ® I n - 
With this substitution, m becomes 

A = -(I + D* ^A 1 ")" 1 [fit "^SU + QX - y 4 ) - (1 <g> I„)TXjv] . (6.7) 

Exactly as in Lemma l4~Tl we have ||fi* _1 || 00 < 2 by (P2), and (I+fit _1 A T ) _1 || 00 < 2 
by (A3). As in the analysis of the state dynamics, all the terms on the right side of 
UEZD are bounded by c||y||oo for a suitable choice of c. Hence, we have ||A||oo ^ c||y||oo- 
Now consider the 6th component of VT*[X, U, A] = y, which can be written 

SjXi + R;Ui + B^Ai = ygj, 1 < z < A. 


Previously, we have shown the existence of a constant c, independent of A, such that 
Moo < c||y ||go and HA^oo < cllyl^, 1 < i < A. By (A2) the smallest eigenvalue 
of Ri is bounded from below by a. Hence, there also exists a constant c such that 
||U i||oo < c||y||oo. Finally, the 3rd component of VT*[X, U, A] = y can be written 

N 

A 0 - TXjv - ^w^ATAi + QiX, + SiUO = y 3 . 

i—l 

By the uniform bounds on HXjUoo, ||U 4 ||oo, and HAjJloo, 1 < i < A, there also exists 
a constant c, independent of A, such that ||A 0 ||oo < c||y||oo- W e summarize these 
results as follows: 

Lemma 6.1. If (A2)-(A3) and (P1)-(P4) hold, then there exists a constant c, 
independent of N, such that the solution of VT*[X, U, A] = y satisfies 

||X||oo + IIUIU + ||A||oo < cllylloo. 


The proof of Theorem 11.11 follows exactly as in [T6]. By Lemma 16.11 p = 
||VT(X*,U*, A*) —1 Hoo is bounded uniformly in A. Choose e small enough that ep < 
1. When we compute the difference VT(X, U, A) - VT(X*,U*,A*) for (X, U, A) 
near (X*,U*, A*) in the oo-norm, the D and D+ constant terms cancel, and we are 
left with terms involving the difference of derivatives of f or C up to second order at 
nearby points. By assumption, these second derivatives are uniformly continuous on 
the closure of Q and on a ball around x*(l). Hence, for r sufficiently small, we have 


||VT(X,U,A)- VT(X*,U*,A*)|| 00 <e 
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whenever 


maxIHX-XlcoJU-iriUJA- A*|U} < r. (6.8) 

By Lemma ETT1 it follows that ||T(X*,U*, A*)|| < (1 — ^e)r//r for all N sufficiently 
large. Hence, by Proposition ^. II there exists a solution to T(X, U, A) = 0 satisfying 
(EH). Moreover, by m and (EH, the estimate (11.121) holds. To complete the proof, 
we need to show that (X, U) is a local minimizer for (EB- After replacing the KKT 
multipliers by the transformed quantities given by m , the Hessian of the Lagrangian 
is a block diagonal matrix with the following matrices forming the diagonal blocks: 

Wi ^(i,u) ff (Xi>Ui, Aj), 1 <i<N, 

WiVf x ‘ u) H(X i , U*, AO + V^ a . j1l) C'(X i ), i = IV, 

where H is the Hamiltonian. In computing the Hessian, we assume that the X and U 
variables are arranged in the following order: Xi, Ui, X 2 , U 2 , ..., Xjv, Ujv. By (A2) 
the Hessian is positive definite when evaluated at (X*, U*, A*). By continuity of the 
second derivative of C and f, and by the convergence result (11.121) . we conclude that 
the Hessian of the Lagrangian, evaluated at the solution of T(X, U, A) = 0 satisfying 
(16.81) . is positive definite for N sufficiently large. Hence, by the second-order sufficient 
optimality condition ES Thm. 12.6], (X, U) is a strict local minimizer of (12.11) . This 
completes the proof of Theorem 11.11 

7. Numerical illustration. Let us consider the unconstrained control problem 
previously introduced in m- 

min {— x(2) : x(t) = |(— x{t) + x[t)u(t) — u(t) 2 ), x(O) = 1} (7.1) 

The optimal solution and associated costate are 

x*(t) = 4/a(f), a(t) = 1 + 3 exp(2.5f), 
u*(t) = x*(t)/ 2, 

A *(t) = — a 2 {t) exp(—2.5t)/[exp(—5) + 9exp(5) + 6]. 

Figure mi plots the logarithm of the sup-norm error in the state, control, and costate 
as a function of the number of collocation points. Since these plots are nearly linear, 
the error behaves like clO~ aN where a « 0.6 for either the state or the control 
and a « 0.8 for the costate. In Theorem 11.11 the dependence of the error on N is 
somewhat complex due to the connection between rj and N. As we increase N , we 
can also increase 77 when the solution is infinitely differentiable, however, the norm 
of the derivatives also enters into the error bound as in (EH- Nonetheless, in cases 
where the solution derivatives can be bounded by c v for some constant c, it is possible 
to deduce an exponential decay rate for the error as observed in [T21 Sect. 2], 

8 . Conclusions. A Radau collocation scheme is analyzed for an unconstrained 
control problem. For a problem with a smooth solution and a Hamiltonian which 
satisfies a strong convexity assumption at a local minimizer of the continuous problem, 
we show that the discrete approximation has a local minimizer in a neighborhood 
of the continuous solution, and as the number of collocation points increases, the 
distance in the sup-norm between the discrete solution and the continuous solution is 
0(N 2 ~ ri ) when the continuous solution has 77 + 1 continuous derivatives, 77 > 3, and 
the number of collocation points N is sufficiently large. A numerical example is given 
which exhibits an exponential convergence rate. 
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Fig. 7.1. The base 10 logarithm of the error in the sup-norm as a function of the number of 
collocation points. 


9. Appendix. Before stating property (P3) in the Introduction, we showed that 
D* is an invertible matrix. In this section, we give an analytic formula for the inverse. 
Proposition 9.1. The inverse o/D* is given by 


D\'■ 1 = ui N Mj(l) + ff Mj(r)dT, 1 < i < N, 1 < j < N, 

Dfj^ 1 = cv N , 1 <i<N, 

D x N -f = UN Mj{ 1), 1 <j<N, 

where Mj, 1 < j < N, is the Lagrange interpolating basis relative to the point set t\, 
.. tjv- 1 - That is, 

N-l 

M i( r )=n I ^r L ’ j = N-l. 

, l~j l~i 
2—1 J 
i^J 


Proof. The relation (11.111) holds for any polynomial p of degree at most N — 1. 
Let p £ denote the vector with i-th component p{ji). In vector form, the system 
of equations (11.111) can be expressed D-^p = p — ejsip(T) /u>n ■ Multiply by D-^ _1 and 
exploit the identity D-f _1 ejv = — oon! of (16.61) to obtain 

_1 P = p- lp(l). (9.1) 

Since p is a polynomial of degree at most N— 2, we can only specify the derivative 
of p at N — 1 distinct points. Given any j satisfying 1 < j < N, let us insert in (19.11) 
a polynomial p £ Vn-i satisfying 

p(rj) = 1 and p(Tj) = 0 for all i < N, i ^ j. 

A specific polynomial with this property is 
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Since p n = p(l) = 0, the last component of the right side of (19.11) vanishes to give the 
relation D]^ 1 + Dy~^p( 1) = 0. In (16.61) we showed that all the elements in the last 
column of D* _1 are equal to — ojn , and by (19.21) . p(l) = 1). Hence, we obtain the 

relation 

^T 1 = -£)^- 1 p(l)=a; i vp(l)=a; J vM j (l), 1 < j < TV. (9.3) 

Finally, let us consider D\~ 1 for i < N and j < TV. We combine the *-th 
component of m for * < TV with (19.21) to obtain 

(D t ~ 1 p) i =J^ Mj{r )dr . (9.4) 

Recall that all components of p vanish except for the j-th, which is 1, and the TV-th, 
which is Mj{ 1) by (19.21) . Hence, (19.41) and the fact that the elements in the last column 
of Dt - 1 are all — u>n yield 

D]- 1 = Mj(T)dT - D\x l Mj{l) = wjvMf(l) + J Mj(r)dr 
This completes the proof. □ 

As noted in the Introduction, the elements in the last row of are positive 

and sum to 2, and the last row of D^W -1 / 2 is positive and has Euclidean norm 
\/2. Numerically, we observe that the absolute row sums for the first TV — 1 rows 
of are always less than 2, while the Euclidean norm of the first TV — 1 rows 

of D( '. N ~W 1 / 2 are always less than \/2. Properties (P3) and (P4) are less tight in 
the sense that the bounds only hold as equalities asymptotically. Tables 19.11 and 19.21 
show ||D* - 1 || 00 and the maximum Euclidean norm of the rows of Dt _1 W _1 / 2 for 
an increasing sequence of dimensions. 


TV 

25 

50 

75 

100 

125 

150 

norm 

1.995376 

1.998844 

1.999486 

1.999711 

1.999815 

1.999871 

TV 

175 

200 

225 

250 

275 

300 

norm 

1.999906 

1.999928 

1.999943 

1.999954 

1.999962 

1.999968 


Table 9.1 
||Dt -iHoo 


TV 

25 

50 

75 

100 

125 

150 

norm 

1.412209 

1.413691 

1.413982 

1.414083 

1.414130 

1.414156 

TV 

175 

200 

225 

250 

275 

300 

norm 

1.414171 

1.414181 

1.414188 

1.414193 

1.414196 

1.414199 


Table 9.2 

Maximum Euclidean norm for the rows of [W 1 / 2 D^] —1 
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